!-------------------------------------------------------------LICENSE--------------------------------------------------------------!
!                                                                                                                                  !
!The MAP code is written in Fortran language for magnetohydrodynamics (MHD) calculation with the adaptive mesh refinement (AMR)    !
!and Message Passing Interface (MPI) parallelization.                                                                              !
!                                                                                                                                  !
!Copyright (C) 2012                                                                                                                !
!Ronglin Jiang                                                                                                                     !
!rljiang@ssc.net.cn                                                                                                                !
!585 Guoshoujing Road. Pudong, Shanghai, P.R.C. 201203                                                                             !
!                                                                                                                                  !
!This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License         !
!as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.             !
!                                                                                                                                  !
!This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of    !
!MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.                        !
!                                                                                                                                  !
!You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software     !
!Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.                                                   !
!                                                                                                                                  !
!-------------------------------------------------------------LICENSE--------------------------------------------------------------!

!==============================================================================================================================================|
recursive subroutine block_generate (current)
!==============================================================================================================================================|

   use parameters
   use block_define

   implicit none

   type(block), pointer, intent(inout) :: current

   integer(4) :: i, i_p, i_block, i_level, i_neigh, sig, i_pp1, i_pm1
   real(8) :: lmt, denominator
   type(block), pointer :: child
!----------------------------------------------------------------------------------------------------------------------------------------------|

!----------------------------------------------------------------------------------------------------------------------------------------------|
!  before regrid to finer mesh, check the level of neighbor of current block.
!  regrid the neighbor block if the level of neighbor is lower the current block.
!----------------------------------------------------------------------------------------------------------------------------------------------|
   do i_neigh = 1, n_neighs
      if (current%neigh_lv(i_neigh) .eq. 0) cycle
      child => current%neigh(i_neigh, current%neigh_lv(i_neigh))%pblock
      if (child%lv .lt. current%lv) call block_generate (child)
   enddo

   allocate (current%child(ratio))
   current%nc = ratio

   n_blocks = n_blocks + current%nc

!----------------------------------------------------------------------------------------------------------------------------------------------|
!     interpolation
!----------------------------------------------------------------------------------------------------------------------------------------------|
   do i_block = 1, current%nc
      child => current%child(i_block)

      write (child%id, "('B', i9.9)") id_index
      child%lv = current%lv + 1
      child%nx = nx_block
      child%p_pos = i_block
      child%dx = current%dx / 2.0d0
      child%dt = current%dt / 2.0d0
      child%t = current%t
      child%pp => current

      allocate (child%x(child%nx), child%ro(child%nx), child%mx(child%nx), child%my(child%nx), child%mz(child%nx),                 &
         child%bx(child%nx), child%by(child%nx), child%bz(child%nx), child%en(child%nx),                                           &
         child%flux_bnd(n_neighs, sys), child%neigh(n_neighs, n_levels))

      id_index = id_index + 1

      do i = 1, child%nx
         i_p = (ng + i + 1) / 2 + (child%p_pos - 1) * nx_no_bnd / ratio
         sig = (-1) ** (i + ng)
         i_pp1 = i_p + 1
         i_pm1 = i_p - 1

         lmt = child%dx / 2.0d0
         child%x(i) = current%x(i_p) + sig * lmt

         denominator = current%ro(i_pp1) - current%ro(i_pm1)
         if (denominator .eq. 0) denominator = eps
         lmt = max((current%ro(i_pp1) - current%ro(i_p)) * (current%ro(i_p) - current%ro(i_pm1)), 0.0d0) / denominator / 2.0d0
         child%ro(i) = current%ro(i_p) + sig * lmt

         denominator = current%mx(i_pp1) - current%mx(i_pm1)
         if (denominator .eq. 0) denominator = eps
         lmt = max((current%mx(i_pp1) - current%mx(i_p)) * (current%mx(i_p) - current%mx(i_pm1)), 0.0d0) / denominator / 2.0d0
         child%mx(i) = current%mx(i_p) + sig * lmt

         denominator = current%my(i_pp1) - current%my(i_pm1)
         if (denominator .eq. 0) denominator = eps
         lmt = max((current%my(i_pp1) - current%my(i_p)) * (current%my(i_p) - current%my(i_pm1)), 0.0d0) / denominator / 2.0d0
         child%my(i) = current%my(i_p) + sig * lmt

         denominator = current%mz(i_pp1) - current%mz(i_pm1)
         if (denominator .eq. 0) denominator = eps
         lmt = max((current%mz(i_pp1) - current%mz(i_p)) * (current%mz(i_p) - current%mz(i_pm1)), 0.0d0) / denominator / 2.0d0
         child%mz(i) = current%mz(i_p) + sig * lmt

         denominator = current%bx(i_pp1) - current%bx(i_pm1)
         if (denominator .eq. 0) denominator = eps
         lmt = max((current%bx(i_pp1) - current%bx(i_p)) * (current%bx(i_p) - current%bx(i_pm1)), 0.0d0) / denominator / 2.0d0
         child%bx(i) = current%bx(i_p) + sig * lmt

         denominator = current%by(i_pp1) - current%by(i_pm1)
         if (denominator .eq. 0) denominator = eps
         lmt = max((current%by(i_pp1) - current%by(i_p)) * (current%by(i_p) - current%by(i_pm1)), 0.0d0) / denominator / 2.0d0
         child%by(i) = current%by(i_p) + sig * lmt

         denominator = current%bz(i_pp1) - current%bz(i_pm1)
         if (denominator .eq. 0) denominator = eps
         lmt = max((current%bz(i_pp1) - current%bz(i_p)) * (current%bz(i_p) - current%bz(i_pm1)), 0.0d0) / denominator / 2.0d0
         child%bz(i) = current%bz(i_p) + sig * lmt

         denominator = current%en(i_pp1) - current%en(i_pm1)
         if (denominator .eq. 0) denominator = eps
         lmt = max((current%en(i_pp1) - current%en(i_p)) * (current%en(i_p) - current%en(i_pm1)), 0.0d0) / denominator / 2.0d0
         child%en(i) = current%en(i_p) + sig * lmt
      enddo

!----------------------------------------------------------------------------------------------------------------------------------------------|
!     neighbors
!----------------------------------------------------------------------------------------------------------------------------------------------|

      do i_neigh = 1, n_neighs
      do i_level = 1, n_levels
         child%neigh(i_neigh, i_level)%pblock => null()
      enddo
      enddo

      if (i_block .eq. 1) then
!----------------------------------------------------------------------------------------------------------------------------------------------|
!     neighbor front
!----------------------------------------------------------------------------------------------------------------------------------------------|
         do i_level = 1, n_levels
            if (.not. associated (current%neigh(1, i_level)%pblock)) cycle
            child%neigh(1, i_level)%pblock => current%neigh(1, i_level)%pblock
            child%neigh_lv(1) = i_level
            current%neigh(1, i_level)%pblock%neigh(2, child%lv)%pblock => child
            current%neigh(1, i_level)%pblock%neigh_lv(2) = child%lv
         enddo
!----------------------------------------------------------------------------------------------------------------------------------------------|
!     neighbor behind
!----------------------------------------------------------------------------------------------------------------------------------------------|
         child%neigh(2, child%lv)%pblock => current%child(i_block + 1)
         child%neigh_lv(2) = child%lv

      elseif (i_block .eq. current%nc) then

!----------------------------------------------------------------------------------------------------------------------------------------------|
!     neighbor front
!----------------------------------------------------------------------------------------------------------------------------------------------|
         child%neigh(1, child%lv)%pblock => current%child(i_block - 1)
         child%neigh_lv(1) = child%lv
!----------------------------------------------------------------------------------------------------------------------------------------------|
!     neighbor behind
!----------------------------------------------------------------------------------------------------------------------------------------------|
         do i_level = 1, n_levels
            if (.not. associated (current%neigh(2, i_level)%pblock)) cycle
            child%neigh(2, i_level)%pblock => current%neigh(2, i_level)%pblock
            child%neigh_lv(2) = i_level
            current%neigh(2, i_level)%pblock%neigh(1, child%lv)%pblock => child
            current%neigh(2, i_level)%pblock%neigh_lv(1) = child%lv
         enddo

      else

         child%neigh(1, child%lv)%pblock => current%child(i_block - 1)
         current%neigh_lv(1) = child%lv
         child%neigh(2, child%lv)%pblock => current%child(i_block + 1)
         current%neigh_lv(2) = child%lv

      endif
   enddo

   return
end subroutine block_generate
